################################################################################
# Script Name: Power Analysis - General Linear Mixed Model Repeated Measures
#
# Purpose: Determine minimum sample size for Generative AI and 
#          Survey Development experiment.
#
# Author: Erica Ann Metheney
#
# Contact: data@gld.gu.se
#
# Notes: 
#        Following Tutorial by Arend & Shäfer 2019, Psychological Methods
#
#
################################################################################

# LOAD PACKAGES
library(simr)

# SET SEED
set.seed(5028)

# RUN LOOPS

#...Range of Values for Simulation Space
alpha.vals = c(0.05,0.05/11)
N.clus.vals = seq(from = 200, to = 350,by = 50)
effect.size.vals.ME = c(0.1,0.41,0.7)
effect.size.vals.IE = c(0.05,0.1,0.2)
results = data.frame(alpha = numeric(),
                     Nclusters = numeric(),
                     EffectSize = numeric(),
                     LowerCI = numeric(),
                     UpperCI = numeric())

#...Model Direct Effect
for(i in 1:length(alpha.vals)){
  for(k in 1:length(N.clus.vals)){
    for(l in 1:length(effect.size.vals.ME)){
      alpha.S <- alpha.vals[i]
      Size.clus <- 6
      N.clus <- N.clus.vals[k]
      
      # CREATE VARIABLES FOR SIMULATION
      g <- as.factor(1:N.clus)
      trt.mat <- data.frame(matrix(c(0,0,0,0,1,0,0,0,1,1,0,0,1,1,0,1,0,1),nrow = 6,byrow = TRUE))
      design.mat = do.call("rbind", replicate(N.clus, trt.mat, simplify = FALSE))
      design.mat$q = rep(g,each = Size.clus)
      
      # ADAPT STANDARDIZED PARAMETERS
      varL1 <- 2
      varL2 <- 1.5
      
      # CREATING CONDITIONAL VARIANCES
      V1 <- varL2
      
      # CREATING A POPULATION MODEL FOR SIMULATION
      b <- c(0,effect.size.vals.ME[l],0.41,0.41,0.1,0.1)
      model <- makeGlmer(y~ X1 + X2 + X3 + X1:X2 + X1:X3 + (1|q),family = binomial(link = "logit"), fixef = b,VarCorr = V1,data = design.mat)
      
      # IMPLEMENT POWER ANALYSIS
      sim.ef <- powerSim(model,fixed("X1"), alpha = alpha.vals[i], nsim = 1000)
      sim.ef.sum <- summary(sim.ef)
      
      # SAVE RESULTS
      results.row = data.frame(alpha = alpha.vals[i],
                               Nclusters = N.clus.vals[k],
                               EffectSize = effect.size.vals.ME[l],
                               LowerCI = sim.ef.sum$lower,
                               UpperCI = sim.ef.sum$upper)
      results = rbind(results,results.row)
    }
  }
}
setwd("C:/Users/xmeter/University of Gothenburg/GLD Projects - Documents/Generative AI and Survey Development/Analyses/Power Analysis")
saveRDS(results,"PowerAnalysisResults_DirectEffectModel.rds")


#...Persona Direct Effect
results = data.frame(alpha = numeric(),
                     Nclusters = numeric(),
                     EffectSize = numeric(),
                     LowerCI = numeric(),
                     UpperCI = numeric())
for(i in 1:length(alpha.vals)){
    for(k in 1:length(N.clus.vals)){
      for(l in 1:length(effect.size.vals.ME)){
        alpha.S <- alpha.vals[i]
        Size.clus <- 6
        N.clus <- N.clus.vals[k]
        
        # CREATE VARIABLES FOR SIMULATION
        g <- as.factor(1:N.clus)
        trt.mat <- data.frame(matrix(c(0,0,0,0,1,0,0,0,1,1,0,0,1,1,0,1,0,1),nrow = 6,byrow = TRUE))
        design.mat = do.call("rbind", replicate(N.clus, trt.mat, simplify = FALSE))
        design.mat$q = rep(g,each = Size.clus)
        
        # ADAPT STANDARDIZED PARAMETERS
        varL1 <- 1
        
        # CREATING CONDITIONAL VARIANCES
        V1 <- varL2
        
        # CREATING A POPULATION MODEL FOR SIMULATION
        b <- c(0,0.41,effect.size.vals.ME[l],0.41,0.1,0.1)
        model <- makeGlmer(y~ X1 + X2 + X3 + X1:X2 + X1:X3 + (1|q),family = binomial(link = "logit"), fixef = b,VarCorr = V1,data = design.mat)
        
        # IMPLEMENT POWER ANALYSIS
        sim.ef <- powerSim(model,fixed("X2"), alpha = alpha.vals[i], nsim = 1000)
        sim.ef.sum <- summary(sim.ef)
        
        # SAVE RESULTS
        results.row = data.frame(alpha = alpha.vals[i],
                                 Nclusters = N.clus.vals[k],
                                 EffectSize = effect.size.vals.ME[l],
                                 LowerCI = sim.ef.sum$lower,
                                 UpperCI = sim.ef.sum$upper)
        results = rbind(results,results.row)
      }
    }
}
setwd("C:/Users/xmeter/University of Gothenburg/GLD Projects - Documents/Generative AI and Survey Development/Analyses/Power Analysis")
saveRDS(results,"PowerAnalysisResults_DirectEffectPersona.rds")


#...Interaction Effect
results = data.frame(alpha = numeric(),
                     Nclusters = numeric(),
                     EffectSize = numeric(),
                     LowerCI = numeric(),
                     UpperCI = numeric())
for(i in 1:length(alpha.vals)){
    for(k in 1:length(N.clus.vals)){
      for(l in 1:length(effect.size.vals.IE)){
        alpha.S <- alpha.vals[i]
        Size.clus <- 6
        N.clus <- N.clus.vals[k]
        
        # CREATE VARIABLES FOR SIMULATION
        g <- as.factor(1:N.clus)
        trt.mat <- data.frame(matrix(c(0,0,0,0,1,0,0,0,1,1,0,0,1,1,0,1,0,1),nrow = 6,byrow = TRUE))
        design.mat = do.call("rbind", replicate(N.clus, trt.mat, simplify = FALSE))
        design.mat$q = rep(g,each = Size.clus)
        
        # ADAPT STANDARDIZED PARAMETERS
        varL1 <- 1
        
        # CREATING CONDITIONAL VARIANCES
        V1 <- varL2
        
        # CREATING A POPULATION MODEL FOR SIMULATION
        b <- c(0,0.41,0.41,0.41,effect.size.vals.IE[l],0.1)
        model <- makeGlmer(y~ X1 + X2 + X3 + X1:X2 + X1:X3 + (1|q),family = binomial(link = "logit"), fixef = b,VarCorr = V1,data = design.mat)
        
        # IMPLEMENT POWER ANALYSIS
        sim.ef <- powerSim(model,fixed("X1:X2"), alpha = alpha.vals[i], nsim = 1000)
        sim.ef.sum <- summary(sim.ef)
        
        # SAVE RESULTS
        results.row = data.frame(alpha = alpha.vals[i],
                                 Nclusters = N.clus.vals[k],
                                 EffectSize = effect.size.vals.IE[l],
                                 LowerCI = sim.ef.sum$lower,
                                 UpperCI = sim.ef.sum$upper)
        results = rbind(results,results.row)
      }
    }
}
setwd("C:/Users/xmeter/University of Gothenburg/GLD Projects - Documents/Generative AI and Survey Development/Analyses/Power Analysis")
saveRDS(results,"PowerAnalysisResults_InteractionEffect.rds")